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Abstract 


In  this  report  we  describe  a  method  for  modelling  the  scattering  from  three- 
dimensional  bathymetric  features  on  the  seabed.  In  particular,  the  features 
are  taken  to  be  azimuthally  symmetric  about  the  z-axis.  This  approach  could 
model  seamount  or  basin  type  features  or  at  higher  frequencies,  smaller  pro¬ 
tuberances  and  depressions.  The  method  uses  a  set  of  point  sources  within 
the  scattering  region  and  exterior  to  the  region  to  satisfy  the  appropriate  con¬ 
tinuity  and  boundary  conditions  along  scattering  features  and  seabed/water 
interface.  Multi-frequency  and  Fourier-synthesized  time  domain  responses  for 
the  scattering  from  such  features  are  presented.  In  addition  we  present  a  semi- 
analytic  modal  solution  which  can  be  used  as  a  benchmark  for  a  bathymetric 
feature  in  a  halfspace  with  only  a  change  in  density 


Resume 


Dans  le  present  rapport,  nous  decrivons  une  methode  pour  la  modelisation  de 
la  diffusion  a  partir  de  caracteristiques  bathvmetriques  tridimensionnelles  du 
plancher  oceanique.  Plus  particulierement,  ces  caracteristiques  sont  clioisies 
de  maniere  a  ce  qu’elles  soient  symetriques  azimutalement  par  rapport  a  l’axe 
des  “z”.  Cette  methode  permet  la  modelisation  de  monts  sous-marins  ou  de 
bassins  ou,  a  des  frequences  superieures,  de  protuberances  ou  de  depressions. 
Cette  methode  utilise  un  ensemble  de  sources  ponctuelles  situees  a  l’interieur 
et  a  l’exterieur  de  la  zone  de  diffusion  afin  de  respecter  les  conditions  adequates 
de  continuite  et  de  limite  le  long  de  caracteristiques  de  diffusion  et  a  l’interface 
plancher  oceanique-eau.  Des  reponses  adimension  temporelle,  multifrequences, 
par  syntliese  de  Fourier,  pour  la  diffusion  par  de  telles  caracteristiques  elements, 
sont  presentees.  Nous  presentons  egalement  une  solution  module  semi-analytique 
qui  peut  etre  utilisee  comnie  reference  pour  une  earaeterist  icjue  batliymetrique 
dans  un  demi-espace  ou  il  n’y  a  qu’une  seule  variation  de  densite. 
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Executive  summary 


BACKGROUND 

The  modelling  and  understanding  of  acoustic  scattering  from  bathymetric  fea¬ 
tures  is  important  in  a  number  of  sonar  applications  where  the  scattering 
of  acoustic  energy  from  seabed  features  produces  reverberation  and  highlight 
features  which  may  interfere  with  the  detection  of  objects  of  interest.  These 
objects  may  be  submarines  or  torpedoes  within  the  water  column  or  sea  mines 
on  the  bottom  or  buried  in  the  underlying  sediment.  Bathymetric  features 
can  also  cast  acoustic  shadows  affecting  the  detection  of  objects  on  the  seabed 
and  they  can  significantly  alter  the  transmission  of  sonar  energy  into  the  sed¬ 
iment.  This  report  presents  a  method  for  computing  the  acoustic  scattering 
from  simple  three-dimensional  bathymetric  features  on  the  seabed. 

SIGNIFICANCE  OF  RESULTS 

It  is  shown  that  the  scattering  from  three-dimensional  bathymetric  features 
(with  an  azimuthal  symmetry)  can  be  modeled  accurately  and  efficiently  for 
low  to  moderate  frequencies  (e.g.,  a  few  kHz).  The  model  allows  one  to  accu¬ 
rately  predict  the  effects  of  simple  bathymetric  features  on  the  scattering  and 
propagation  of  a  sonar  signal.  The  three-dimensional  solutions  which  can  be 
obtained  from  the  methods  of  this  report  should  prove  useful  as  benchmark 
solutions  for  other  scattering  models. 

FUTURE  WORK 

In  the  future,  more  general  and  efficient  methods  (although,  perhaps  more 
approximate)  which  could  be  used  for  more  general  seabed  bathymetry  and 
at  higher  frequencies  than  considered  in  this  report,  should  be  investigated. 
The  methods  described  in  this  report  could  provide  accurate  and  challenging 
benchmark  cases. 

Fawcett,  J.  Modelling  scattering  from  simple  three-dimensional  bathymetric 
features  using  wavefield  superposition.  DRDC  Atlantic  TM  2006-086,  Defence 
R&D  Canada  -  Atlantic. 
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CONTEXTE 

La  modelisation  et  la  comprehension  de  la  diffusion  acoustique  a  partir  de 
caracteristiques  bathymetriques  sont  importantes  pour  un  grand  nombre 
duplications  aux  sonars,  dans  lesquelles  la  diffusion  de  l’energie  acoustique 
par  des  caracteristiques  du  plancher  oceanique  produit  de  la  reverberation  et 
fait  ressortir  des  caracteristiques  qui  peuvent  interferer  lors  de  la  detection  des 
objets  recherches,  Ces  objets  peuvent  etre  des  sous-marins  ou  des  torpilles 
se  trouvant  dans  la  colonne  d’eau  ou  des  mines  marines  deposees  au  fond  ou 
enfouies  dans  les  sediments  sous-jacents.  Les  caracteristiques  bathymetriques 
peuvent  egalement  produire  des  ombres  acoustiques  qui  influent  sur  la  detection 
d’objets  situes  sur  le  plancher  oceanique  et  elles  peuvent  modifier 
considerablement  1a,  transmission  de  l’impulsion  emise  jusque  dans  les  sediments 
Dans  le  present  rapport,  on  presente  une  methode  de  traitement  informatique 
de  la  diffusion  acoustique  produite  par  des  caracteristiques  bathymetriques 
tridimensionnelles  simples  du  plancher  oceanique. 

PORTEE 

II  est  demontre  que  la  diffusion  par  des  caracteristiques  bathymetriques  tridi¬ 
mensionnelles  (a  symetrie  azimutale)  peut  etre  modelisee  avec  precision  et 
efficacement  pour  des  frequences  de  faibles  a  moderees  (de  quelques  kHz).  Le 
modele  permet  de  predire  avec  precision  les  effets  de  caracteristiques 
bathymetriques  simples  sur  la  diffusion  et  la  propagation  du  signal  d’un  sonar. 
Les  solutions  tridimensionnelles  pouvant  etre  obtenues  par  les  metliodes  figu¬ 
rant  dans  le  present  rapport  devraient  etre  utiles  comnie  solutions  de  reference 
pour  d’autres  modeles  de  diffusion. 

RECHERCHES  FUTURES 

Dans  l’avenir,  des  metliodes  plus  generates  et  efficaces  (bien  que,  peut-etre 
plus  approximates )  pouvant  etre  utilisees  pour  la  bathymetrie  plus  general e 
et  pour  des  frequences  plus  elevees  que  celles  considerees  dans  le  present  rap¬ 
port  devraient  etre  etudiees.  Les  metliodes  decrites  dans  le  present  rapport 
pourraient  servir  de  references  precises  et  rigoureuses. 
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1  INTRODUCTION 


The  modelling  of  acoustic  scattering  from  range-dependent  bathymetry  or  sur¬ 
faces  is  a  subject  of  much  interest  in  underwater  acoustics.  For  example, 
bathymetric  features  can  cause  significant  sonar  returns,  reverberation,  which 
can  interfere  with  the  detection  of  a  target  (e.g.,  submarine,  torpedo,  mine). 
A  bathymetric  feature  can  cast  an  acoustic  shadow  on  the  seabed  which  may 
hinder  the  detection  of  mines  located  within  this  shadow.  In  addition,  when 
searching  for  a  buried  target,  the  sonar  signal  which  is  transmitted  into  the 
underlying  sediment  may  be  significantly  altered  by  the  seabed  features.  At 
higher  sonar  frequencies,  these  effects  can  be  caused  by  relatively  small-sized 
features.  Thus  the  ability  to  predict  the  effects  of  seabed  bathymetry  on  sonar 
returns  is  important  to  understanding  sonar  system  performance  in  a  number 
of  applications. 

There  have  been  many  methods  proposed  for  modelling  propagation  and  scat¬ 
tering  in  three-dimensional,  possibly  range-dependent  media.  These  include 
Parabolic  Equation  methods  [1],  coupled  mode  methods  [2-4],  boundary  el¬ 
ement  methods  [5],  multipole  expansion  or  wavefield  superposition  methods 
[6-8]  and  finite  difference  methods  [9].  These  various  methods  all  have  their 
advantages  and  disadvantages  in  terms  of  accuracy,  computation  efficiency  and 
generality.  In  this  report,  we  consider  a  scattering  chamber /wavefield  super¬ 
position  method  [10].  The  scattering  chamber  refers  to  the  fact  that  a  small 
region  surrounding  the  scattering  feature  of  interest  is  defined  and  this  region 
is  coupled  to  the  surrounding  medium  by  a  set  of  interior  and  exterior  virtual 
sources.  We  will  consider  three-dimensional  bathymetric  features  which  are 
azimuthally-symmetric. 

The  advantages  of  the  method  of  this  report  is  that  it  reduces  the  compu¬ 
tational  domain  to  a  region  which  encloses  the  scattering  surface  while  still 
insuring  that  the  resulting  scattered  fields  satisfy  the  boundary  conditions 
of  the  waveguide  outside  the  scattering  region.  This  type  of  approach  was 
used  in  Ref.10  to  model  scattering  from  a  rough  sea  surface.  Wave  super¬ 
position  methods  have  also  been  used  to  model  target  scattering  in  a  wave- 
guide[7].  Schmidt[8]  has  combined  the  wave  superposition  method,  Finite 
Element  modelling  and  the  OASES  propagation  code  to  obtain  a  computa¬ 
tional  method  for  partially-buried  targets  in  a  waveguide.  Our  approach  for 
bathymetric  scattering  is  conceptually  similar  to  the  target  scattering  case.  We 
introduce  an  additional  curve  to  make  the  bathymetric  feature  into  a  closed 
region  which  has  the  properties  of  one  of  the  surrounding  media.  We  repre¬ 
sent  the  field  interior  to  this  region  with  a  set  of  point  sources  (rather  than 
using  a  finite  element  model)  to  characterize  the  region’s  interior  fields  and 
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we  use  another  set  of  point  sources  in  the  interior  of  the  region  to  represent 
the  exterior  field.  These  point  sources  in  the  interior  are  used  in  conjunction 
with  the  Green’s  Function  for  the  surrounding  medium,  thus  insuring  that  the 
appropriate  boundary  conditions  (e.g.,  along  the  water/seabed  interface)  are 
satisfied.  The  point  source  coefficients  are  determined  from  a  system  of  equa¬ 
tions  expressing  the  continuity  conditions  across  the  boundary  defining  the 
bathymetry  region.  Once  the  point  source  coefficients  have  been  determined, 
the  pressure  field  can  be  computed  for  an  arbitrary  receiver  position.  In  addi¬ 
tion,  we  will  consider  azimuthallv-symmetric  features.  This  allows  us  to  solve 
the  full  three-dimensional  scattering  problem  as  a  sequence  of  two-dimensional 
problems. 
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2  THEORY 


For  the  pressure  field,  p,  and  the  Green’s  Function  G  in  an  azimuthallv  sym¬ 
metric  waveguide  we  can  write, 


p{r\  f§  9) 
Pinc(r,  z,  9) 
G{r,  z,  9 ;  r',  z' ,  9’) 


J ~2pm{r,z )  cos  (m3) 

m= 0 
oo 

PmC(r,z)  cos (m,9) 

m= 0 
oo 

71  ffm(r,  Z, ;  r\  z’)  cos {m{9  -9')), 


(1) 


m= 0 


where  (r,  ,2)  is  a  two-dimensional  cylindrical  coordinate  system  and  9  is  the  az¬ 
imuthal  angle  about  the  axis  of  symmetry,  the  z-axis.  The  terms  gm(r,  Z, ;  r',  z') 
are  computed  from  the  azimuthal  transform  of  the  Green’s  Function.  In  this 
report  we  will  be  considering  both  the  Green’s  Function  for  three-dimensional 
free  space  and  the  Green’s  Function  for  a  half  space  and  their  azimuthal  trans¬ 
forms.  For  example,  the  half-space  Green’s  Function  for  the  mth  azimuthal 
order  can  be  written  for  the  source  point  ( r',z ')  in  the  upper  half  space, 


9m(r,z,\rftz')  =  /  Jm.{pr)Jm{pr') 


x  (exp(j7ijf  ~  ^|)  +  R{p)  e^Pl +  z')))pdp  Z>Q 

hi 

f°°  T  ,  \  T  t  /xr(p)exp(-?;(72^  +  7i^))  , 

=  /  Jm{pr)Jm{pr) - : - pdp,z<  0 

Jo  Hi 


(2) 


where  71  =  \J 'h1  / cf  —  p2  and  R(p)  and  T(p)  are  the  reflection  and  transmission 
coefficients  for  the  incident  field  in  the  upper  halfspace.  A  similar  expression 
exists  for  the  case  that  the  source  points  are  in  the  lower  half  space  with  7  ,R(p) 
and  T(p)  replaced  by  the  appropriate  expressions. 


The  Green’s  Function  in  a  free  space  is  given  by 


G(r,  z,  6;  r\  9')  = 


exp(icoR) 

R 


where 

R  =  yj r 2  +  r'2  —  2 rr'  cos(0  —  6')  +  (z  —  z')2. 


(3) 

(4) 
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[Many  authors  include  a  factor  of  47T  in  the  denominator  of  Eq.(3),  but  to 
simply  our  notation  and  to  be  consistent  with  Eq.(2),  we  will  not  include  this 
factor].  The  azimuthal  transform  of  this  function  can  be  expressed  as 


gm{r,z,;r%z')  = 


exp  (i'jjR) 

R 


cos  (mij’)  dil’ 


(5) 


where  ip  =  6  —  91 ,  recalling  from  Eq.(4)  that  R  depends  on  9  and  9'  through 
their  difference.  There  are  a  variety  of  analytical  expressions  for  Eq.(5)  either 
as  a  series  of  Legendre  polynomials  and  spherical  Hankel  or  Bessel  functions 
or  in  terms  of  elliptical  integrals  [11].  Also,  one  can  efficiently  compute  the 
expression  for  several  values  of  m  using  FFT  methods.  The  tradeoff,  in  this 
case,  is  the  requirement  to  store  the  resulting  values  for  subsequent  computa¬ 
tion.  For  this  report,  we  simply  discretize  the  integral  of  Eq.(5)  and  perform 
the  integration  numerically  as  required. 

For  the  halfspace  Green’s  Function,  we  can  rewrite  it  in  the  form,  for  z'  >  0 


gm{r,z,\r',z!)  = 


x 


Jm{pr)Jm{pr') 


( R{P )  ~  RgS)ex y{  'ni{z  +  z'))) 

ih 


pdp  +  tdir  +  tref ,  ^  >  0 


(6) 


'Hi 


Tas  exp (  — 772(3  -  z')) 
172 

tf.rans  1  <  0 


]  pdp 


(7) 


where 


tdir  — 
tref  —  Ras 


^  exp [ii/jrjr2  +  r 12  —  2 rr'  cos  ip  +  (z  —  z')2) 
o  \J  r2  +  r'2  —  2  rr1  cos  V’  +  (z  —  z')2 

H  exp {iujyjr2  +  r'2  —  2 rr1  cos  ip  +  (z  +  z')2) 


t 


trans 


—  T 

—  -L  n  . 


'o  \Jr2  +  r'2  —  2 rr'  cos  (/’)  +  ( z  +  z')2 
exp (iu>2\/r2  +  r 12  —  ’2rr'  Los  v  +  (z  —  z')2) 
o  \Jr2  +  r'2  —  2 rr'  cos  A  +  (z  —  z')2 


cos  (mp)d'i 
cos  (m'ip)dip 


cos  (imp)  dp’  (8) 
(9) 


and 


Ra.s  — 


T  — 

±  as  — 


p2  ~  Pi 
p2  +  Pi 
2P2 

P2  +  Pi 


(10) 
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The  terms  of  Eq.(9)  represent  singular  terms  which  occur  for  z'  — >■  z,  r'  — >■  r 
and  for  z  =  0  (for  trej  and  ttrans)-  In  the  case  of  a  Boundary  Integral  Equa¬ 
tion  Method  (BIEM),  the  singular  behaviour  of  t(k,r  .trej  and  ttrans  would  need 
to  be  analytically  accounted  for.  For  the  method  of  images  considered  in 
this  paper,  there  will  always  be  a  non-zero  distance  between  the  source  and 
the  surface  where  the  boundary  conditions  are  applied.  Similar  expressions 
to  those  of  Eqs.(9)  and  (10)  are  obtained  for  z'  <  0.  In  Fig.  1  we  show  a 
schematic  cross  section  of  a  “bump”  on  the  seafloor.  In  this  case  the  feature 
is  pointing  into  the  upper  medium.  The  solid  line  is  the  actual  bathymetry 
and  the  dashed  line  represents  an  arbitrary  numerical  boundary  placed  into 
the  bottom.  The  squares  represent  sources  placed  a  specified  distance  ex¬ 
terior  to  the  region  and  are  used  to  represent  the  field  within  the  region 
Q.  These  are  free  space  sources  (suitable  transformed  into  the  azimuthal  do¬ 
main)  with  the  sound  speed  c2.  The  density  within  the  regions  is  that  of  the 
lower  halfspace  p2.  The  triangles  represent  sources  shifted  slightly  to  the  in¬ 
terior  of  the  bounding  curve  C  and  are  associated  with  the  half-space  Green’s 
Functions  (in  the  azimuthally-transformed  domain)  for  these  particular  source 
positions.  The  resultant  exterior  field  satisfies  the  boundary  conditions  along 
the  water/seabed  interface.  In  the  case  that  we  are  considering  a  basin  (i.e.,  a 
negative  deformation  of  the  seabed)  the  situation  is  almost  exactly  the  same, 
except  that  now  the  interior  values  for  the  sound  speed  and  density  for  Q  are 
those  of  the  upper  half  space.  If  a  deformation  is  both  positive  and  negative 
with  respect  to  the  flat  surface  z  =  0m,  then  there  will  be  an  area  of  the 
exterior  regions  which  is  not  the  simple  half-space  (for  example,  an  intrusion 
of  the  upper  medium  into  the  bottom  half  space).  In  this  case,  one  must  use  a 
larger  bounding  surface  and  introduce  additional  unknown  point  sources  along 
the  actual  surface  deformation.  We  will  not  consider  this  case  in  this  paper. 


Taking  N  sources  with  unknown  amplitudes  {a,;}  in  the  exterior  region  (repre¬ 
senting  the  interior  field)  and  N  sources  in  the  interior  region  bj. (representing 
the  exterior  field),  the  continuity  equations  across  the  bounding  curve  become, 
first,  for  the  continuity  of  the  pressure  field, 

N  N 

Y  a*GLA,  Xi  +  s,)  =  Y  hGZ(xk- xr  -6i)+  pZc(xk)  (11) 

i-  I  i=  1 

and  the  normal  derivative  across  the  boundary 


1  y-  dG{n(Xk\  Xi  +  Si) 

Pin  4?  fl?;  dn 

i=i 


Pext(Xk)  ,=1 


(yw££»W;V  -5)  ,  acc(W) 


dn 


dn 


) 


(12) 

In  these  equations  we  have  indicated  the  discrete  points  along  the  curve  C  as 
X i  =  ( /’, .  2, )  and  Xk  =  {rk,Zk).  The  source  points  are  displaced  slightly  from 
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Figure  1:  A  schematic  of  a  positive  bathymetric  feature.  The  scattering  chamber 
region  Q  is  surrounded  by  the  2  curves  Y  |  which  is  the  surface  of  the  feature  and 
r2  which  is  an  additional  curve  introduced  into  the  seabed  to  enclose  Q. 
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these  points,  along  the  normal,  in  the  positive  direction  +5i  for  the  exterior 
sources  and  —  d,  for  the  interior  sources.  We  have  used  the  superscript  / 
with  the  Green’s  Function  in  Eqs.(ll)  and  (12)  to  denote  free  space  and  the 
superscript  w  to  denote  the  surrounding  waveguide.  In  this  report  we  consider 
a  2  halfspace  waveguide  but  the  same  method  can  also  be  used  for  the  case 
that  surrounding  waveguide  is  bounded  above  by  a  free  surface. 

In  the  case  that  the  bathymetric  feature  is  a  positive  deformation  into  the 
water  column,  then  pin  =  p2:  for  a  negative  deformation  pin  =  p\.  The 
density  term  pext(Xk)  is  equal  to  pi  for  such  that  Zk  >  0  and  is  equal 
to  p2  otherwise.  It  should  also  be  noted  that  although  we  are  stressing  the 
modeling  of  scattering  by  three-dimensional  bathymetric  features,  the  region 
can  also  model  an  embedded  region  of  different  sound  velocity  and  density. 
In  this  case  the  interior  values  of  c  and  p  would  be  different  from  those  of  the 
2  surrounding  halfspaces. 
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3  A  SEMI-ANALYTICAL 
BENCHMARK  CASE 


In  previous  references  [12,13]  modal  solutions  were  derived  for  acoustic  half- 
space  problems  in  the  case  that  the  2  halfspaces  had  the  same  sound  speeds  but 
different  densities.  In  Ref.  13  this  allowed  for  the  solution  of  two-dimensional 
scattering  from  an  embedded  cylinder  in  terms  of  a  coupled  mode  solution.  In 
this  paper,  we  will  derive  a  corresponding  three-dimensional  solution.  In  par¬ 
ticular,  we  will  describe  the  coupled  mode  problem  to  solve  for  each  azimuthal 
order  m.  Let  us  consider  a  sphere  which  is  exactly  half-buried  between  two 
half  spaces,  each  with  sound  speed  cext  and  densities  pi  and  p2 .  The  sphere 
has  an  internal  sound  speed  c,;„  and  density  pjn.  In  the  case  that  c,;n  =  cext 
and  p.in  =  pi  we  have  a  basin  problem  and  for  pin  =  p2  we  have  a  positive 
deformation  problem.  The  general  solution  in  the  interior  of  the  sphere  can 
be  expanded  in  terms  of  the  angular  functions, 

iC(4>,R)  =  KPZ(m s«.  (13) 

where  P™  denotes  a  Legendre  Function,  and  N™  is  the  normalization  constant 
such  that 

r  1 

/  — (if’™(cos  4>))2  sin  4>cl(j)  =  1.  (14) 

JO  Pin 

This  angular  function  has  the  associated  radial  function,  jn(kintR )  where  jn 
denotes  a  spherical  Bessel  Function. 

In  the  exterior  domain,  the  angular  functions  have  the  form  for  n  =  0,  2,  4, ..., 


C(C0S  <j>)  =  Pn+m( COS  <j>)  0  <  0  <  7T 


and  for  n  =  1,  3,  5, ..., 

Pn+mi C°S^)  0<  <£<7t/2 

^P^m(C0S^)  7t/2<0<7T. 

pi 

These  modal  functions  are  then  normalized  such  that 


C(cos</>)  = 
C(cos  <f>)  = 


(15) 


(16) 


P(^) 


7"_ 


')2(cos  <p)  sin  odo  =  1. 


(17) 


The  associated  radial  functions  for  these  angular  functions  are  hn(kextR).  The 
density  term  p(cp)  in  Eq.(17)  is  simply  pi  for  0  <  <p  <  7t/2  and  p2  for  tt/2  < 
0  <  7r.  The  specified  form  of  the  angular  modal  functions  is  used  because 
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for  n  +  m  =  rn,  m  +  2,  m  +  4, ...  the  associated  Legendre  functions  are  even 
(and  have  zero  derivative  with  respect  to  <j>)  and  thus  are  continuous  across 
f:  =  0,  (f>  ~  7r/2  and  trivially  satisfy 


1  dr  1  dr 

pi  d(f)  p2  dcf) 


4>  =  7t/2. 


(18) 


The  Legendre  functions  P,™m(cos0)  are  odd  with  respect  to  <p.  at  <p  =  tt/2  for 
n  =  1,  3,  5, ...  and  thus,  in  fact,  the  definition  of  r  as  in  Eq.(16)  is  continuous 
(zero)  at  <p  =  tt/2  and  satisfies  the  boundary  condition  of  Eq.(18).  We  can 
now  set  up  a  system  of  equations  to  be  satisfied  along  the  surface  (in  this 
coordinate  system,  R  =  a  and  0  <  <f>  <  tt)  of  the  sphere, 


N  N 

y  Cii'i (cos  <!>)ji+m{kinR)  =  y  d-iTi (cos  (j))hi+m{kextR) 

i= 0  i= 0 

N 

T  ^  ^  2ikextTi{cos  (f>)ji-i.m(kextR)hj/+m(kextRs)Ti(cos  <f>s)  (19) 

i= o 


1 

Pin 


N 

y  cd’ii cos  (p)kinjl+m  {kinR) 

i= 0 


1 

Pext ( 0 ) 


N 

^  djTj(cos  (p) keXfh,-^m[keXfR) 

i= 0 


N 

T  ^  ^  2zTj(cos  <p)kex^j,j^m[kextR)hi-^m[keXfRs)Ti\cos  <ps)^j.  (20) 

j=0 


For  the  system  of  equations  in  Eq.(19),  we  project  l/pin'0j(cos  </>)  sin  0,  j  = 
0, ...,  JV  on  both  sides  and  integrate.  For  the  systems  of  equations  in  Eq.(20)  we 
project  Tj  (cos  cp)  sin  d>  on  both  sides  and  integrate.  The  last  terms  in  Eqs.(19) 
and  (20)  are  the  expansion  of  the  incident  field  from  a  point  source  (of  the 
form  exp (ikR)/R)  located  (in  spherical  coordinates)  at  (Rs,(ps,9  =  0).  The 
resulting  system  of  equations  has  the  matrix/vector  form 


(  D1  R, 

\  R>2  D2 


(21) 


Here  D\  and  D2  are  diagonal  matrices  with  entries  ji+m(kinR )  and  hi+m(kextR ) 
respectively.  The  matrices  R\  and  /?2  have  the  form 

1  r 

Ri(p,q)  =  hq+m(kextR) —  /  r/)(coso!r,,(coso)siiir'jr/r'j.  (22) 

Pin  J  0 

1  r 

R2(p,q )  =  jq+m(kinR) —  /  ( '(/(cos  o)  r;/(cos  o)  sill  ado.  (23) 

Pin  J  0 
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The  solution  of  the  system  of  equations  of  Eq.(21)  yields  the  interior  coeffi¬ 
cients  c.  and  the  exterior  coefficients  cl  which  can  then  be  used  to  compute  the 
scattered  psc(r,z\'m)  and  total  pressure  fields  p(r,z;m).  The  scattered  pres¬ 
sure  field  is  defined  as  the  difference  between  the  total  field  and  the  incident 
field  which  would  exist  in  the  absence  of  the  boundary  feature.  These  two- 
dimensional  fields  are  computed  for  M  values  of  the  azimuthal  order  m  and 
the  three-dimensional  fields  are  then  computed  (for  example,  the  scattered 
field)  from 

M 

Psc(r ,  z,  0)  =  ^  empsc{r,  z\  m )  cos (m0)  (24) 

m= 0 

where  em  =  1  for  m  =  0  and  em  =  2  otherwise. 
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4  NUMERICAL  COMPUTATIONS 


The  system  of  equations  for  the  more  general  case  of  the  two  halfspaces  having 
different  sound  velocities,  Eqs.(ll)  and  (12)  is  very  similar  to  Eqs.(19)  and  (20) 
but  now  we  use  the  conditions  at  the  N  values  of  Xk  to  solve  for  the  coefficients 
{a,;}  and  {&*}  of  the  displaced  point  sources.  Once  these  coefficients  have 
been  computed,  they  can  be  used  in  conjunction  with  the  appropriate  Green’s 
Function  to  compute  the  scattered  and  total  fields  within  the  2  halfspaces. 
This  is  done  for  each  azimuthal  order  m.  The  final  solution  is  computed  from 
the  azimuthal  sums  of  the  form  of  Eq.(l).  In  cases  where  the  incident  pressure 
field  is  added  to  the  solution  (for  example,  to  compute  the  total  field  in  the 
exterior  region),  the  incident  field  is  computed  only  once  in  the  form, 

poo 

Pmc{r ,  z;  rs,  zz)  =  /  J0{p\r  -  rs \ )K{p;  zs)pdp  (25) 

J  o 

where  K(p;zs )  represents  the  appropriate  kernel  for  the  upper  or  lower  half¬ 
space.  In  other  words,  we  compute  the  incident  field  with  a  three-dimensional 
expression  with  the  origin  at  the  source  position. 

For  the  benchmark,  density-jump  modal  case,  the  scattering  region  was  a 
half-buried  sphere.  One  could  allow  for  different  fractions  of  burial  by  using- 
functions  based  upon  a  displaced  centre  of  the  sphere.  However,  in  this  case 
the  coupling  equations  would  be  more  complicated  and  the  surface  of  the 
region  would  not  be  a  constant  range  surface  for  this  new  centre.  Thus,  for 
the  benchmark  modal  solution  we  will  only  consider  the  half-buried  situation. 
However,  the  Scattering  Chamber  Wavefield  Superposition  (SOWS)  approach 
can  be  used  for  more  arbitrarily-shaped  surfaces. 

4.1  Density-jump  benchmark  solution 

As  a  first  numerical  benchmark  we  will  use  both  the  benchmark  modal  ap¬ 
proach  and  the  SCWS  approach.  We  will  consider  a  spherical  region  (radius 
=  1.0  m)  with  an  interior  sound  speed  cp  =  1500  m/s  and  p  =  1  g/cm3.  The 
upper  surrounding  half  space  has  a  sound  speed  of  1500  m/s  and  a  density  of 
1  g/cm 3  and  the  lower  half  space  a  sound  speed  of  1500  m/s  and  a  density 
of  5  g/cm3.  This  scenario  represents  a  half-spherical  basin  intruding  into  the 
lower,  high-density  half  space.  A  point  source  is  located  at  a  range  of  10  m  and 
an  angle  of  5°  from  the  z-axis.  The  bounding  chamber  is  simply  a  half-circle 
in  the  (r,  ,?)-space  in  this  case.  A  total  of  122  points  were  used  to  discretize 
this  curve.  The  external  and  internal  point  sources  (in  the  azimuthal  domain) 
were  displaced  a  distance  of  0.05  m  along  the  normal  direction.  The  resulting 
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system  of  equations  to  solve  is  244  x  244.  The  number  of  azimuthal  terms 
used  was  given  by  the  expression  nm  =  27t//1500  x  1.5  +  9.  As  a  benchmark 
solution  we  use  the  modal  solution,  using  the  same  number  of  azimuthal  terms 
as  for  the  scattering-chamber  approach.  The  number  of  modes  considered  for 
each  azimuthal  order  is  given  by  the  expression,  nmofje  =  27t//1500  x  3  +  16.  In 
Fig.  2  we  show  a  comparison  of  the  computed  total  fields  along  the  line  z  =  0 
(solid  line  is  modal  solution,  dashed  is  scattering  chamber  approach).  As  can 
be  seen,  the  agreement  is  excellent.  The  points  x  =  ±1  are  the  boundaries  of 
the  basin.  In  Fig.  3  we  show  the  same  computations  at  2500  Hz.  Once  again, 
it  can  be  seen  that  the  agreement  with  the  benchmark  solution  is  excellent. 
In  Fig.  4  we  show  a  grid  of  pressure  field  values  computed  using  the  scattering 
chamber  approach.  As  can  be  seen,  this  “basin”  has  a  significant  effect  upon 
the  propagation  and  there  is  a  strong  focusing  feature  near  the  bottom.  It  is 
also  important  to  note  that  there  are  no  observable  discontinuities  in  going 
across  the  boundary  r  =  1  where  r  denotes  the  distance  from  the  origin.  The 
values  for  x  <  0  are  computed  by  using  6  =  tt  in  the  azimuthal  sum. 

The  pressure  field  is  predicted  to  be  somewhat  ill-behaved  near  x  =  ±1  m ,  Z  = 
0  m.  This  is  because  there  is  a  discontinuity  in  the  normal  (r)  derivative  at 
this  point  because  of  the  discontinuous  boundary  conditions  which  arise  on  the 
boundary  of  the  scattering  chamber  across  the  boundary  z  =  0.  In  Fig.  5  we 
show  the  real  part  of  the  amplitudes  for  the  exterior  (the  first  122  points)  and 
interior  sources  (the  last  122  points)  for  m  =  1  and  /  =  2500 Hz.  As  can  be 
seen,  the  solutions  become  larger  and  more  erratic  towards  :  0.  However, 

as  seen  from  Fig.  4  they  combine  to  give  a  nicely  continuous  solution  across 
the  boundary  r  =  1. 

For  the  modal  method,  the  natural  coordinate  system  is  a  spherical  system  and 
the  pressure  grid  computations  (using  MATLAB)  were  most  easily  vectorizable 
in  such  a  coordinate  system.  For  the  SCWS  approach  the  best  coordinate 
system  is  a  (r,z)  coordinate  system.  In  this  approach,  the  computation  of 
the  halfspace  Green’s  Function  using  Eqs.(7)  and  (9)  explicitly  defined  the 
asymptotic  terms  involved  with  the  density  jump.  The  remainder  wavenumber 
integrals  are  zero  in  this  case  of  only  a  density  jump.  However,  in  order  to 
simplify  the  vectorizing  of  the  grid  computations  once  the  coefficients  have 
been  determined,  the  wavenumber  integral  approach  (for  example,  Eq.(7)  with 
the  determined  coefficients  for  the  external  field)  is  used  with  the  asymptotic 
and  direct  terms  included.  Some  care  must  be  taken  with  these  integrals  (in 
terms  of  the  upper  limit  of  the  integral  and  its  offset  into  the  complex-p  plane) 
to  ensure  that  meaningful  results  are  obtained  for  very  short  source/receiver 
ranges  in  the  computations. 
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x(m) 

Figure  2:  A  comparison  of  the  computed  pressure  files  for  the  basin  scattering 
feature  -  500  Hz,  modal  solution  (solid  line)  and  scattering  chamber  method 
(dashed  line) 


x(m) 

Figure  3:  A  comparison  of  the  computed  pressure  files  for  the  basin  scattering 
feature  -  2500  Hz,  modal  solution  (solid  line)  and  scattering  chamber  method 
(dashed  line) 
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Figure  4:  Computed  two-dimensional  vertical  slice  of  the  three-dimensional 
pressure  field  using  the  scattering  chamber  approach  at  2500  Hz.  The  basin 
feature  is  evident  and  the  incident  field  is  due  to  a  point  source  located  10  m 
away  and  at  5°  off  vertical 
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Figure  5:  The  real  part  of  the  solution  for  m  =  1  for  the  point  source 
amplitudes  for  the  scattering  chamber  approach  at  2500  Hz.  The  values  near  61 
and  183  represent  the  solution  near  r  =  1  m  and  z  =  0. 


4.2  A  general  seabed  feature 

We  now  consider  the  more  general  situation  of  an  axisymmetric  feature  with 
the  surrounding  lialfspace  consisting  of  one  fluid  c\,pi  overlying  the  seabed 
C-2 , p2  ■  The  method  we  propose  is  also  easily  extendable  to  the  case  of  waveguide 
where,  for  example,  there  might  be  an  upper  bounding  pressure  release  surface. 
In  this  case,  we  simply  replace  the  half  space  Green’s  functions  for  the  exterior 
region  with  the  appropriate  waveguide  Green’s  function.  We  consider  the 
following  analytical  expression  for  the  deformation  of  the  seabed, 

.  ,  gw  , 2-7T r,, 

z(r)  =  -U  +  cos(— ))•  (26) 

Z  uj 

We  will  consider  a  =  2  and  w  =  4;  this  is  a  “bump”  which  has  a  height  of 
2  m  and  a  radial  extent  of  2  m.  It  is  important  to  recall  that  actual  feature 
is  this  feature  rotated  about  the  z-axis,  so  that  the  true  lateral  extent  in  a 
horizontal  slice  is  4  m.  A  point  source  in  the  water  column  at  10  m  range  and 
at  a  height  of  1  m  above  the  seabed  is  considered.  For  the  upper  medium, 
ci  =  1500  m/s ,  p\  =  1  g/cm3  and  for  the  lower  sediment  C2  =  1800  m/s  and 
p2  =  1.5 g/cm3.  The  lower  medium  has  an  attenuation  of  0.3  dB/A.  The 
two-dimensional  slice  of  the  scattering  feature,  the  additional  bounding  curve 
in  the  bottom,  and  the  source  points  are  shown  in  Fig.  6.  The  number  of 
azimuthal  terms  used  in  the  computation  varied  with  respect  to  frequency 
according  to  the  relation  nm  =  67r//1500  +  6.  This  results  in  37  (38,  including 
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m  =  0)  azimuthal  terms  for  the  frequency  of  2500  Hz.  In  the  MATLAB 
matrix  inversion  we  set  the  singular  value  tolerance  at  10-8  as  the  system 
matrix  tended  to  have  conditioning  issues  for  high  azimuthal  orders.  As  with 
the  density  contrast  benchmark  solutions,  there  are  122  interior  and  exterior 
point  sources.  The  distance  of  these  sources  from  the  curve  is  important  in 
obtaining  accurate  solutions  and  it  was  found  that  optimal  position  varied 
somewhat  with  frequency.  For  lower  frequencies,  the  distance  was  greater 
than  at  higher  frequencies.  This  is  reasonable  if  one  considers  the  longer 
wavelength  at  lower  frequencies.  In  the  following  we  compute  the  pressure 
fields  for  frequencies  between  500  and  2500  Hz  in  steps  of  25  Hz.  After  some 
numerical  experimentation  with  these  distances,  we  used  the  expression 

d{f)  =  .12  -  (/  -  500)/750  x  .04  500  Hz  <=  f  <=  1250  Hz 

d(f)  =  0.08  1250  Hz  <  f  <=  2500  Hz  with  (27) 

d(  f)  =  0.07  1725  Hz  <=  f  <=  1925  Hz ,  (28) 

for  our  computations.  This  sequence  of  distances  yielded  good  scattering  re¬ 
sults  as  a  function  of  frequency  over  the  band.  In  fact,  the  first  two  set  of 
distances  in  Eq.(28)  yielded  sufficiently  good  results.  However,  an  inspection 
of  the  computed  scattered  fields  as  a  function  of  frequency  indicated  a  slightly 
spurious  response  at  1725  Hz  in  terms  of  the  backscattered  energy.  Although 
this  anomaly  is  small  and  is  insignificant  in  subsequent  time  domain  compu¬ 
tations,  it  is  an  indication  that  the  computational  parameters  are  not  optimal 
in  this  frequency  regime  and  we  found  that  decreasing  the  distance  improved 
the  results  for  this  frequency  resulting  in  the  point  source  offsets  of  the  third 
case  of  Eq.(28).  This  illustrates  a  problem  with  the  method  for  multifrequency 
computations.  There  are  a  number  of  computational  parameters  which  need 
to  be  specified  for  the  method:  the  number  of  discrete  points,  the  additional 
contour  added  in  either  the  seabed  or  water  column,  the  offsets  of  the  points 
sources  from  the  curves  C\  and  (A,  and  even  the  parameters  used  in  the  nu¬ 
merical  computation  of  the  Green’s  Functions.  It  is  desirable  to  have  simple 
formulae  to  specify  these  parameters  as  the  frequency  is  varied  in  the  numeri¬ 
cal  code  to  automatically  compute  the  pressure  field  as  a  function  of  frequency. 
However,  the  computational  results  did  seem  somewhat  sensitive  to  the  spec¬ 
ification  of  the  point  source  offsets  and  the  parameters  for  the  wavenumber 
integral  for  the  Green’s  functions.  Some  care  and  numerical  experimentation 
was  required  to  determine  good  parameter  values  as  a  function  of  frequency. 

In  Figs.  7  and  8  we  show  the  computed  fields  for  6  different  frequencies  for  a 
vertical  grid  of  receiver  locations  and  for  a  horizontal  grid  at  a  fixed  vertical 
distance  of  z  =  0.55  m.  In  these  figures  there  are  no  observable  discontinuities 
of  the  field  across  the  surface  and  the  additional  bounding  curve  in  the  bot¬ 
tom.  It  is  interesting  to  note  that  since  the  feature  is  penetrable,  some  energy 
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Figure  6:  Two-dimensional  cross-section  of  bathymetric  feature.  The  portion 
above  z  =  0  is  the  actual  feature.  A  smooth  bounding  curve  is  added  in  the 
bottom.  The  blue  plus  signs  are  the  discrete  points  along  the  curves  where  the 
continuity  equations  are  imposed.  The  green  plus  signs  are  the  source  points 
interior  to  the  domain  which  represent  the  scattered  fields  in  the  exterior;  the 
half  space  Green's  Function  (in  the  azimuthal  domain)  are  used  for  these  sources. 
The  maroon  plus  signs  represent  the  fields  within  the  domain  and  these  use  the 
free  space  Green's  Function  for  the  interior  values  of  sound  speed  (complex  in 
the  case  of  attenuation). 
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Figure  7:  Two-dimensional  vertical  (y  =  0  rn)  cross-section  of  the  total  pressure 
field  for  6  different  frequencies:  (a)  500  Hz  (b)  900  Hz  (c)  1300  Hz  (d)  1700  Hz 
(e)  2100  Hz  (f)  2500  Hz.  The  point  source  is  to  the  right  of  the  figure  (off  the 
grid).  The  absolute  values  of  the  levels  are  shown  in  the  interval  [0  0.3]) 


propagates  into  the  bump  and  out  the  other  side.  In  Fig.  8  the  circle  denotes 
the  boundary  of  the  bathymetric  feature  at  this  altitude  above  the  seabed  and 
hence  the  acoustic  parameters  within  this  circle  correspond  to  those  of  the 
seabed.  As  one  would  expect,  the  interference  patterns  between  the  direct 
and  reflected  incident  fields  and  the  scattered  field  become  more  complicated 
for  higher  frequencies. 

In  Fig.  9  we  show  the  variation  of  the  absolute  value  of  the  total  and  scat¬ 
tered  fields  as  a  function  of  frequency  for  a  fixed  point  in  the  grid  (x=1.23  m, 
z=2.05  m)  which  is  a  point  of  moderate  scattering.  As  can  be  seen  the  2  curves 
are  well  behaved  as  a  function  of  frequency.  We  did  find  that  the  scattered 
field  amplitudes  did  show  more  significant  fluctuations  as  a  function  of  fre¬ 
quency,  particularly  for  the  higher  frequencies,  at  fixed  grid  locations  close  to 
the  surface  of  the  feature.  In  this  case,  some  of  these  fluctuations  may  be  due 
to  the  closeness  of  the  receiver  points  to  either  the  interior  or  exterior  point 
sources  combined  with  any  numerical  inaccuracies  in  the  point  source  ampli- 
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Figure  8:  Two-dimensional  horizontal  cross-section  (for  fixed  z  =  0.55  rn)  of  the 
total  pressure  field  for  6  different  frequencies:  (a)  500  Hz  (b)  900  Hz  (c)  1300 
Hz  (d)  1700  Hz  (e)  2100  Hz  (f)  2500  Hz.  The  point  source  is  to  the  right  of  the 
figure  (off  the  grid).  The  absolute  values  of  the  levels  are  shown  in  the  interval 
[0  0.3],  The  circle  indicates  the  boundary  of  the  bathymetric  feature  at 
z  =  0.55  m.) 
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Figure  9:  The  variation  of  the  absolute  value  of  the  total  field  (top  panel)  and 
scattered  field  (bottom  panel)  for  the  point  ( x,y,z )  =  (1.23m,  0.0m,  2.05m) 

tudes  and  numerical  inaccuracies  in  computing  the  resulting  pressure  values 
(particularly  for  locations  very  close  to  the  point  sources). 

In  Figs.  10  -  12  we  show  some  time  domain  computations.  The  frequency 
results  were  transformed  from  the  frequency  domain  to  the  time  domain  by 
using  the  expression, 


N 

p(x,y,z,t )  =  ^{y ^p(x,y,z;  fn )  exp(— ?;2?r fnt)S(fn)}  (29) 

n=l 

where  for  the  source  weighting  we  used 

S(f)  =  exp (  — (/  -  1500)2/(2  x  5002)).  (30) 
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Figure  10:  Two-dimensional  vertical  (y  =  0  m)  cross-section  of  the  total 
pressure  field  for  6  different  times:  (a)  4.5  ms  (b)  5.1  ms  (c)  5.7  ms  (d)  6.3  ms 
(e)  6.9  ms  (f)  7.5  ms.  The  point  source  is  to  the  right  of  the  figure  (off  the 
grid).  The  values  are  shown  in  the  interval  [-5.  5.]) 


In  Fig.  10  the  resulting  pulses  are  shown  for  a  vertical  plane  for  6  instances 
of  time.  The  pulse  is  propagating  from  the  right  to  the  left  of  the  plot.  The 
incident  pulse,  in  the  seabed,  is  slightly  ahead  of  the  pulse  in  the  water  because 
of  the  higher  sound  speed  in  the  seabed.  The  pulse  can  be  seen  to  travel  up 
the  scattering  features  and  also  through  the  feature.  A  reflected  wavefront 
begins  to  appear  in  panel(d). 

In  Fig.  11  the  propagation  of  the  pulse  in  a  horizontal  plane  can  be  seen.  In 
Panel  (c)  the  wavefront  is  just  entering  into  the  feature  at  this  depth  and 
begins  to  pull  ahead  of  the  wavefront  on  the  outside.  There  is  a  significant 
focusing  of  energy  within  the  feature.  In  Fig.  12  we  consider  this  interior 
scattering  with  a  finer  time  sampling  where  a  reflected  wavefront  can  be  seen 
quite  noticeably  within  the  feature.  This  wavefront  exits  the  feature  into  the 
water  but  is  quite  weak  at  this  point. 
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Figure  11:  Two-dimensional  horizontal  (z  =  0.55  rn)  cross-section  of  the  total 
pressure  field  for  6  different  times:  (a)  4.5  ms  (b)  5.1  ms  (c)  5.7  ms  (d)  6.3  ms 
(e)  6.9  ms  (f)  7.5  ms.  The  point  source  is  to  the  right  of  the  figure  (off  the 
grid).  The  values  are  shown  in  the  interval  [-5.  5.].  The  circle  indicates  the 
boundary  of  the  bathymetric  feature  at  z  =  0.55  m.) 
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Figure  12:  Two-dimensional  horizontal  (z  =  0.55  rnj  cross-section  of  the  total 
pressure  field  for  6  different  times:  (a)  7.50  ms  (b)  7.65  ms  (c)  7.80  ms  (d)  7.95 
ms  (e)  8.10  ms  (f)  8.25  ms.  The  point  source  is  to  the  right  of  the  figure  (off 
the  grid).  The  values  are  shown  in  the  interval  [-5.  5.].  The  circle  indicates  the 
boundary  of  the  bathymetric  feature  at  z  =  0.55  m.) 
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5  SUMMARY 


In  this  memorandum  we  have  described  an  approach  for  computing  the  pres¬ 
sure  field  scattered  from  three-dimensional  seabed  features.  The  features  were 
taken  to  be  azimuthally-symmetric  and  the  deformation  to  the  surround¬ 
ing  flat  seabed  was  taken  to  be  of  one  polarity  (positive  -  a  “bump”  and 
negative  a  “basin”).  The  computational  method  solved  a  sequence  of  two- 
dimensional  problems  and  then  used  Fourier  synthesis  to  construct  the  full 
three-dimensional  solution.  A  benchmark  coupled  mode  solution  for  the  case 
of  a  hemispherical  deformation  with  a  density-contrast  seabed  was  considered. 
This  type  of  solution  was  used  to  show  that  the  wave  superposition/scattering 
chamber  approach  was  very  accurate  for  this  case.  This  more  general  method 
was  then  used  to  solve  the  scattering  from  a  2  m  high,  penetrable  “bump”.  The 
pressure  field  was  accurately  computed  for  a  large  range  of  frequencies  for  both 
vertical  and  horizontal  planes.  At  the  highest  frequency,  37  azimuthal  orders 
were  computed.  The  multi-frequency  computations  allowed  for  the  compu¬ 
tation  of  time  domain  solutions  for  an  incident  pulse.  Computationally,  the 
method  is  straightforward.  However,  it  was  found  that  some  care  had  to  be 
taken  with  the  positioning  of  the  point  sources. 
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